####
Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given data set, and should not be used in the context of making policy decisions without external consultation from scientific experts.
According to a recent report, overall tobacco use increased in youths (middle school and high school students) in the United States in 2017 and 2018, despite previous years of declining use.
This major increase is attributed to an increase in the use of electronic cigarette (e-cigarette) products.
E-cigarettes are referred to by many different names, including but not limited to:
The devices vary greatly:
See this CDC guide and the American Lung Association website for more information.
The report found that:
During 2017–2018, current use of any tobacco product increased 38.3% (from 19.6% to 27.1%) among high school students and 28.6% (from 5.6% to 7.2%) among middle school students; e-cigarette use increased 77.8% (from 11.7% to 20.8%) among high school students and 48.5% (from 3.3% to 4.9%) among middle school students.
In 2018, the Federal Drug Administration (FDA) in the United States stated that e-cigarette usage use among youth reached:
“nothing short of an epidemic proportion of growth”
In this case study, we will be investigating the same data used in the report that generated the above findings. This data comes from the The National Youth Tobacco Survey (NYTS).
Gentzke, Andrea S., Melisa Creamer, Karen A. Cullen, Bridget K. Ambrose, Gordon Willis, Ahmed Jamal, and Brian A. King. “Vital Signs: Tobacco Product Use Among Middle and High School Students - United States, 2011-2018.” MMWR. Morbidity and Mortality Weekly Report 68 (6): 157–64 (2019).
Our main question:
What vaping brands and flavors appear to be used the most frequently?
We will base this on the following survey questions:
> “During the past 30 days, what brand of e-cigarettes did you usually use?”
>“What flavors of tobacco products have you used in the past 30 days?”
Have vaping rates possibly influenced tobacco/nicotine use?
In this case study, we will cover how to import data from multiple files efficiently, how to import data from excel files, and how to make a variety of visualizations to compare multiple groups across time. We will also demonstrate how to work with codebooks. We will cover the concept of survey weighting and introduce the srvyr package. We will discuss the difference between pooled cross-sectional data and panel data. We will especially focus on using packages and functions from the Tidyverse for wrangling data, such as tidyr and dpyr and for visualization, such as as ggplot2. The tidyverse is a library of packages created by RStudio. While some students may be familiar with previous R programming packages, these packages make data science in R especially efficient.
We will begin by loading the packages that we will need:
library(here)
library(readxl)
library(stringr)
library(purrr)
library(dplyr)
library(readr)
library(tidyr)
library(ggplot2)
library(forcats)
library(srvyr)
library(cowplot)| Package | Use |
|---|---|
| here | to easily load and save data |
| readxl | to import the data in the excel files |
| stringr | to manipulate the character strings within the data |
| purrr | to import the data in all the different excel and csv files efficiently |
| dplyr | to arrange/filter/select/compare specific subsets of the data |
| readr | to import the CSV file data |
| tidyr | to rearrange data in wide and long formats |
| ggplot2 | to make visualizations with multiple layers |
| forcats | to allow for reordering of factors in plots |
| syrvr | to use survey weights |
| cowplot | to allow plots to be combined |
The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.
According to the cited Morbidity and Mortality Weekly Report this was what was already known about this topic and the implications of this study:
Importantly, the vapors used in e-cigarettes contain harmful chemicals:
There are some important considerations regarding this data analysis to keep in mind:
According to Wikipedia: >Panel data differs from pooled cross-sectional data across time, because it deals with the observations on the same subjects in different times whereas the latter observes different subjects in different time periods
While gender and sex are not actually binary, the data used in this analysis only contains information for groups of individuals who answered the survey questions as male or female.
The data in this case study comes from the National Youth Tobacco Survey (NYTS) which is an annual survey that asks students in high school and middle school (grades 6-12) about tobacco usage in the United States of America.
The data for this survey is freely available online at this website with data from 1999, 2000, 2002, 2004, 2006, 2009, and 2011-2019. We will be using data from 2015-2019 due to the fact that these years are the most recent that asked questions regarding e-cigarette usage.
Each year includes documentation, such as a codebook and an excel file containing the data:
Therefore, since we are using data from 2015-2019, the data we are interested in is located in 5 different excel files, each with their own codebook.
The codebook contains information describing the data within the excel file.
As you can see the excel file contains very short variable names and values, and it is not clear what they mean without the codebook:
The codebook explains what the variables (the columns) are: 
And the codebook explains what the values for each variable are:
We will explain more later about what the values on the right indicate.
The reason that there are codebooks for each year is because the questions asked year varied slightly.
The data in this survey is what is called pooled cross-sectional data. In other words, different subsets of students are surveyed each year and it is not clear which, if any , individuals participate from one year to the next.
Since these excel files are so large (each has roughly 20,000 rows), it takes a bit of time for the data to load. To make the process faster, we previously imported these files, selected only our questions of interest, and saved this data as csv files.
Click here for details on how the data was originally imported
First we created a list of file names of all the different excel files. Using the here() function of the here package, we looked in all the directories of the project. The list.files() function looked for all files with .xlsx within these sub-directories.
[1] "docs/2015-nyts-dataset-and-codebook-microsoft-excel/nyts2015.xlsx"
[2] "docs/2016-nyts-dataset-and-codebook-microsoft-excel/nyts2016.xlsx"
[3] "docs/2017-nyts-dataset-and-codebook-microsoft-excel/nyts2017.xlsx"
[4] "docs/2018-nyts-dataset-and-codebook-microsoft-excel/nyts2018.xlsx"
[5] "docs/2019-nyts-dataset-and-codebook-microsoft-excel/nyts2019.xlsx"All the files were read using read_excel() of the readxl package. Using the map() function of the purrr package this was done efficiently for all of the excel files in the list using one command. The . is used to indicate that we want to apply the read_excel() function to the data that we just piped into the map() function.
Here will also used the %>% pipe which can be used to define the input for later sequential steps.
This created a single list of tibbles (one for each file).
Each excel file name was extracted using the str_extract() function of the stringr package. Here we are keeping occurrences of the character string “nyts201” followed by a “5”,“6”,“7”,“8”, or “9”.
[1] "nyts2015" "nyts2016" "nyts2017" "nyts2018" "nyts2019"These names became the names of the tibbles in the list of tibbles.
Specific columns were selected from each of the tibbles using the variable name, as identified in the codebook for being of interest. In some cases functions like starts_with() of the dplyr package were used to select several variables. Most of the survey questions about tobacco use start with an "E" or a "C" according to the codebooks.
tbl_files[["nyts2015"]] <- tbl_files[["nyts2015"]] %>%
dplyr::select(psu, #Primary Sampling Unit
finwgt,#Analysis Weight
stratum,#Sampling stratum
Qn1, #Age
Qn2, #Sex
Qn3, #Grade
starts_with("E",
ignore.case = FALSE),
starts_with("C",
ignore.case = FALSE),
)
tbl_files[["nyts2016"]] <- tbl_files[["nyts2016"]] %>%
dplyr::select(psu,
finwgt,
stratum,
Q1, #Age
Q2, #Sex
Q3, #Grade
starts_with("E",
ignore.case = FALSE),
starts_with("C",
ignore.case = FALSE),
Q50A, #Menthol # What flavors of tobacco products have you used in the past 30 days? (Select one or more)
Q50B, #Clove or spice
Q50C, #Fruit
Q50D, #Chocolate
Q50E, #Alcoholic Drink
Q50F, #Candy/Desserts/Other Sweets
Q50G, #Some Other Flavor Not Listed Here
Q50H #I Did Not Use Flavored Tobacco Products In the Past
)
tbl_files[["nyts2017"]] <- tbl_files[["nyts2017"]] %>%
dplyr::select(psu,
finwgt,
stratum,
Q1, #Age
Q2, #Sex
Q3, #Grade
starts_with("E",
ignore.case = FALSE),
starts_with("C",
ignore.case = FALSE),
Q50A, #Menthol # What flavors of tobacco products have you used in the past 30 days? (Select one or more)
Q50B, #Clove or spice
Q50C, #Fruit
Q50D, #Chocolate
Q50E, #Alcoholic Drink
Q50F, #Candy/Desserts/Other Sweets
Q50G, #Some Other Flavor Not Listed Here
Q50H #I Did Not Use Flavored Tobacco Products In the Past
)
tbl_files[["nyts2018"]] <- tbl_files[["nyts2018"]] %>%
dplyr::select(psu,
finwgt,
stratum,
Q1, #Age
Q2, #Sex
Q3, #Grade
starts_with("E",
ignore.case = FALSE),
starts_with("C",
ignore.case = FALSE),
Q50A, #Menthol # What flavors of tobacco products have you used in the past 30 days? (Select one or more)
Q50B, #Clove or spice
Q50C, #Fruit
Q50D, #Chocolate
Q50E, #Alcoholic Drink
Q50F, #Candy/Desserts/Other Sweets
Q50G, #Some Other Flavor Not Listed Here
Q50H #I Did Not Use Flavored Tobacco Products In the Past
)
tbl_files[["nyts2019"]] <- tbl_files[["nyts2019"]] %>%
dplyr::select(psu,
finwgt,
stratum,
Q1, #Age
Q2, #Sex
Q3, #Grade
starts_with("E",
ignore.case = FALSE),
starts_with("C",
ignore.case = FALSE),
Q40, #Brand, e-cigarettes
Q62A, #Menthol # What flavors of tobacco products have you used in the past 30 days? (Select one or more)
Q62B, #Clove or spice
Q62C, #Fruit
Q62D, #Chocolate
Q62E, #Alcoholic Drink
Q62F, #Candy/Desserts/Other Sweets
Q62G, #Some Other Flavor Not Listed Here
)A directory was created using the base dir.create() function called data_reduced for the csv files. New csv files were created for each of the tbls in the list using the write_csv() function of the readr package. This was done all at once using the base mappy() function.
#
dir.create("docs/data_reduced")
mapply(write_csv, tbl_files, path=here(paste0("docs/data_reduced/",
names(tbl_files),
'.csv'))
)
#trying to use map...
tbl_files %>%
map(~write_csv(. , path=paste0("docs/data_reduced2/nyts", pull(year)[[1]], ",csv")))
tbl_files %>%
map(~write_csv(. , path=paste0("docs/data_reduced2/nyts", names(.), ",csv")))Now we will show how to read in the data from the five csv files that were created from the five different excel files.
So let’s take a look at our data. We can get a good sense of it using the glimpse() function of the dplyr package.
Rows: 17,711
Columns: 29
$ psu <chr> "015438", "015438", "015438", "015438", "015438", "015438"…
$ finwgt <dbl> 216.7268, 324.9620, 324.9620, 397.1552, 264.8745, 264.8745…
$ stratum <chr> "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "B…
$ Qn1 <dbl> 10, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,…
$ Qn2 <dbl> 2, 1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2…
$ Qn3 <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5…
$ ECIGT <dbl> 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 2, 1, 2, 2, 1, 2, 1, 2, 1…
$ ECIGAR <dbl> 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 2…
$ ESLT <dbl> 2, 2, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ EELCIGT <dbl> 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 1…
$ EROLLCIGTS <dbl> 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2…
$ EFLAVCIGTS <dbl> 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ EBIDIS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ EFLAVCIGAR <dbl> 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1, 1, 2…
$ EHOOKAH <dbl> 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ EPIPE <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ ESNUS <dbl> 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ EDISSOLV <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CCIGT <dbl> 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CCIGAR <dbl> 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CSLT <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CELCIGT <dbl> 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CROLLCIGTS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CFLAVCIGTS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CBIDIS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CHOOKAH <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CPIPE <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CSNUS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
$ CDISSOLV <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…Currently, it isn’t very clear what most of the variables indicate.
We want to rename variables like Qn1 to something more meaningful like Age.
To do this we will use the rename() function of the dplyr package. The new name is always listed first before the =.
We also want to update the values for age and grade, so that they are more understandable.
Recall from the codebook for this year, that age isn’t listed in the way one might expect:
The same is true for grade:
This is why it is so important to always check the codebook!!
We also need to add new variables about usage of brands and specific flavors of e-cigarettes, because although later years have questions about this, this year unfortunately only has broad questions about flavors. Eventually we want to put the data from each year together, so we need the same variables for each year. Thus we need to add variables that just say NA for the brand and flavor data for this year.
To change the existing variables and create new variables, we will use the mutate function of the dplyr package.
nyts_data[["nyts2015"]] <- nyts_data[["nyts2015"]] %>%
dplyr::mutate(
brand_ecig=NA,
menthol=NA,
clove_spice=NA,
fruit=NA,
chocolate=NA,
alcoholic_drink=NA,
candy_dessert_sweets=NA,
other=NA,
no_use=NA) Now we can see how the data has changed:
Rows: 17,711
Columns: 38
$ psu <chr> "015438", "015438", "015438", "015438", "015438"…
$ finwgt <dbl> 216.7268, 324.9620, 324.9620, 397.1552, 264.8745…
$ stratum <chr> "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "BR3",…
$ Age <dbl> 10, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1…
$ Sex <dbl> 2, 1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2, 1, …
$ Grade <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, …
$ ECIGT <dbl> 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 2, 1, 2, 2, 1, …
$ ECIGAR <dbl> 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ ESLT <dbl> 2, 2, 2, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, …
$ EELCIGT <dbl> 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, …
$ EROLLCIGTS <dbl> 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, …
$ EFLAVCIGTS <dbl> 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ EBIDIS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ EFLAVCIGAR <dbl> 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, …
$ EHOOKAH <dbl> 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, …
$ EPIPE <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ ESNUS <dbl> 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ EDISSOLV <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CCIGT <dbl> 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CCIGAR <dbl> 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CSLT <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CELCIGT <dbl> 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CROLLCIGTS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CFLAVCIGTS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CBIDIS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CHOOKAH <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CPIPE <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CSNUS <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ CDISSOLV <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ brand_ecig <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ menthol <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ clove_spice <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ fruit <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ chocolate <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ alcoholic_drink <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ candy_dessert_sweets <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ other <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ no_use <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …For the next years, we do have flavor questions, but not a brand question. We can replace these vague question names like Q50A with values like menthol. The next 3 years have the same structure for these questions. For each of these years, we also want create a brand variable that is just NA values. We can do this efficiently if we create a function to do all of these steps at once. Then we can use the map_at() function of the purrr package to apply the function we created (for renaming variables etc.) to the all of the data except for the data from 2015. By using vars() inside of the map_at() function we can specify what tibbles within our nyts_data list that we want to include or exclude.
Update_survey <- function(x) { x %>%
rename(Age=Q1,
Sex=Q2,
Grade=Q3,
menthol=Q50A,
clove_spice=Q50B,
fruit=Q50C,
chocolate=Q50D,
alcoholic_drink=Q50E,
candy_dessert_sweets=Q50F,
other=Q50G,
no_use=Q50H) %>%
mutate(brand_ecig=NA)
}
#few options to apply to the data:
#nyts_data <-nyts_data %>% map_at(vars(nyts2016, nyts2017, nyts2018), Update_survey)
#nyts_data <-nyts_data %>% map_at(c("nyts2016", "nyts2017", "nyts2018"), Update_survey)
nyts_data <-nyts_data %>% map_at(vars(-nyts2015, -nyts2019), Update_survey)The final year has a slightly different data structure - It actually has a brand_ecig variable.
nyts_data[["nyts2019"]] <- nyts_data[["nyts2019"]] %>%
rename(brand_ecig=Q40,
Age=Q1,
Sex=Q2,
Grade=Q3,
menthol=Q62A,
clove_spice=Q62B,
fruit=Q62C,
chocolate=Q62D,
alcoholic_drink=Q62E,
candy_dessert_sweets=Q62F,
other=Q62G) %>%
mutate(no_use="missing")Now let’s take a look at the variables for each of the years:
$nyts2015
[1] "psu" "finwgt" "stratum"
[4] "Age" "Sex" "Grade"
[7] "ECIGT" "ECIGAR" "ESLT"
[10] "EELCIGT" "EROLLCIGTS" "EFLAVCIGTS"
[13] "EBIDIS" "EFLAVCIGAR" "EHOOKAH"
[16] "EPIPE" "ESNUS" "EDISSOLV"
[19] "CCIGT" "CCIGAR" "CSLT"
[22] "CELCIGT" "CROLLCIGTS" "CFLAVCIGTS"
[25] "CBIDIS" "CHOOKAH" "CPIPE"
[28] "CSNUS" "CDISSOLV" "brand_ecig"
[31] "menthol" "clove_spice" "fruit"
[34] "chocolate" "alcoholic_drink" "candy_dessert_sweets"
[37] "other" "no_use"
$nyts2016
[1] "psu" "finwgt" "stratum"
[4] "Age" "Sex" "Grade"
[7] "ECIGT" "ECIGAR" "ESLT"
[10] "EELCIGT" "EHOOKAH" "EROLLCIGTS"
[13] "EFLAVCIGAR" "EPIPE" "ESNUS"
[16] "EDISSOLV" "EBIDIS" "CCIGT"
[19] "CCIGAR" "CSLT" "CELCIGT"
[22] "CHOOKAH" "CROLLCIGTS" "CPIPE"
[25] "CSNUS" "CDISSOLV" "CBIDIS"
[28] "menthol" "clove_spice" "fruit"
[31] "chocolate" "alcoholic_drink" "candy_dessert_sweets"
[34] "other" "no_use" "brand_ecig"
$nyts2017
[1] "psu" "finwgt" "stratum"
[4] "Age" "Sex" "Grade"
[7] "ECIGT" "ECIGAR" "ESLT"
[10] "EELCIGT" "EHOOKAH" "EROLLCIGTS"
[13] "EPIPE" "ESNUS" "EDISSOLV"
[16] "EBIDIS" "CCIGT" "CCIGAR"
[19] "CSLT" "CELCIGT" "CHOOKAH"
[22] "CROLLCIGTS" "CPIPE" "CSNUS"
[25] "CDISSOLV" "CBIDIS" "menthol"
[28] "clove_spice" "fruit" "chocolate"
[31] "alcoholic_drink" "candy_dessert_sweets" "other"
[34] "no_use" "brand_ecig"
$nyts2018
[1] "psu" "finwgt" "stratum"
[4] "Age" "Sex" "Grade"
[7] "ECIGT" "ECIGAR" "ESLT"
[10] "EELCIGT" "EHOOKAH" "EROLLCIGTS"
[13] "EPIPE" "ESNUS" "EDISSOLV"
[16] "EBIDIS" "CCIGT" "CCIGAR"
[19] "CSLT" "CELCIGT" "CHOOKAH"
[22] "CROLLCIGTS" "CPIPE" "CSNUS"
[25] "CDISSOLV" "CBIDIS" "menthol"
[28] "clove_spice" "fruit" "chocolate"
[31] "alcoholic_drink" "candy_dessert_sweets" "other"
[34] "no_use" "brand_ecig"
$nyts2019
[1] "psu" "finwgt" "stratum"
[4] "Age" "Sex" "Grade"
[7] "ECIGT" "ECIGAR" "ESLT"
[10] "EELCIGT" "EHOOKAH" "EROLLCIGTS"
[13] "EPIPE" "ESNUS" "EDISSOLV"
[16] "EBIDIS" "EHTP" "CCIGT"
[19] "CCIGAR" "CSLT" "CELCIGT"
[22] "CHOOKAH" "CROLLCIGTS" "CPIPE"
[25] "CSNUS" "CDISSOLV" "CBIDIS"
[28] "CHTP" "brand_ecig" "menthol"
[31] "clove_spice" "fruit" "chocolate"
[34] "alcoholic_drink" "candy_dessert_sweets" "other"
[37] "no_use" It’s looking better!
We also want to replace the value of 19 for Age to be ">18" and the value of 13 for Grade to be replaced with "Ungraded/Other" Also according to the codebooks, numeric values of 1 indicate a survey answer of FALSE, while a value of 2 indicates TRUE. Sex needs to be recoded, but it is a bit unclear if it is encoded slightly differently across years.
Let’s create a function to make all these updates except for Sex:
Update_values <- function(x) { x %>%
mutate(Age=as.numeric(Age)+8,
Grade=as.numeric(Grade)+5)%>%
mutate(Age=as.factor(Age),
Grade=as.factor(Grade),
) %>%
mutate_all(~ replace(., . %in% c("*", "**"), NA))%>%
mutate(Age=recode(Age,
`19` = ">18",
),
Grade=recode(Grade,
`13` = "Ungraded/Other")) %>%
mutate_at(vars(starts_with("E", ignore.case = FALSE),
starts_with("C", ignore.case = FALSE)),
list(~recode(.,
`1` = TRUE,
`2` = FALSE,
.default = NA,
.missing = NA)))
}
nyts_data <-nyts_data %>% map(.,Update_values)Now let’s update the Sex encoding:
It’s a bit difficult to tell if the encoding was the same across years, so let’s check using the count() function of the dplyr package.
According to the codebook, we should have:
1) 8,958 males in 2015 2) 10,438 males in 2016 3) 8,881 males in 2017
4) 10,069 males in 2018
5) 9,803 males in 2019
$nyts2015
# A tibble: 3 x 2
Sex n
<dbl> <int>
1 1 8958
2 2 8622
3 NA 131
$nyts2016
# A tibble: 3 x 2
Sex n
<dbl> <int>
1 1 10438
2 2 10082
3 NA 155
$nyts2017
# A tibble: 3 x 2
Sex n
<dbl> <int>
1 1 8881
2 2 8815
3 NA 176
$nyts2018
# A tibble: 3 x 2
Sex n
<dbl> <int>
1 1 10069
2 2 9920
3 NA 200
$nyts2019
# A tibble: 3 x 2
Sex n
<chr> <int>
1 .N 116
2 1 9803
3 2 9099Thus, it looks like males were encoded as 1 for each year.
Update_sex <- function(x){x %>%
mutate(Sex = as.factor(Sex)) %>%
mutate(Sex = recode(Sex,
`1`= "male",
`2` = "female"))
}
nyts_data <-nyts_data %>% map(., Update_sex)
count_sex <-function(x){x %>%count(Sex)}
nyts_data %>% map(., count_sex)$nyts2015
# A tibble: 3 x 2
Sex n
<fct> <int>
1 male 8958
2 female 8622
3 <NA> 131
$nyts2016
# A tibble: 3 x 2
Sex n
<fct> <int>
1 male 10438
2 female 10082
3 <NA> 155
$nyts2017
# A tibble: 3 x 2
Sex n
<fct> <int>
1 male 8881
2 female 8815
3 <NA> 176
$nyts2018
# A tibble: 3 x 2
Sex n
<fct> <int>
1 male 10069
2 female 9920
3 <NA> 200
$nyts2019
# A tibble: 3 x 2
Sex n
<fct> <int>
1 .N 116
2 male 9803
3 female 9099Looks good!
The years (2016-2019) that have flavors also need the flavor data to be logical:
Update_flavors <- function(x){x %>%
mutate_at(vars(menthol:no_use),
list(~recode(.,
`1` = TRUE,
.default = FALSE,
.missing = FALSE))) }
nyts_data <- nyts_data %>% map_at(vars(-nyts2015), Update_flavors)Now there are just a few changes needed that are specific to 2019:
nyts_data[["nyts2019"]] <- nyts_data[["nyts2019"]] %>%
mutate_all(~ replace(., . %in% c(".N",".S",".Z"), NA)) %>%
mutate(psu=as.character(psu)) %>%
mutate(brand_ecig = recode(brand_ecig,
`1` = "Other", #levels 1,8 combined to `Other`
`2` = "Blu",
`3` = "JUUL",
`4` = "Logic",
`5` = "MarkTen",
`6` = "NJOY",
`7` = "Vuse",
`8` = "Other"))Great! Now we have all the same variables and our values don’t need to be handled any differently for any of the years. Thus we can combine the data across years.
Rows: 95,465
Columns: 41
$ year <chr> "nyts2015", "nyts2015", "nyts2015", "nyts2015", …
$ psu <chr> "015438", "015438", "015438", "015438", "015438"…
$ finwgt <dbl> 216.7268, 324.9620, 324.9620, 397.1552, 264.8745…
$ stratum <chr> "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "BR3",…
$ Age <fct> 18, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, …
$ Sex <fct> female, male, male, male, female, female, male, …
$ Grade <fct> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, …
$ ECIGT <lgl> FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRU…
$ ECIGAR <lgl> TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FA…
$ ESLT <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EELCIGT <lgl> FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRU…
$ EROLLCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EFLAVCIGTS <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, F…
$ EBIDIS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ EFLAVCIGAR <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, F…
$ EHOOKAH <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ EPIPE <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ ESNUS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EDISSOLV <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CCIGT <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, …
$ CCIGAR <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, …
$ CSLT <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CELCIGT <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, …
$ CROLLCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CFLAVCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CBIDIS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CHOOKAH <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CPIPE <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CSNUS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CDISSOLV <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ brand_ecig <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ menthol <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ clove_spice <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ fruit <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ chocolate <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ alcoholic_drink <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ candy_dessert_sweets <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ other <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ no_use <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ EHTP <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ CHTP <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …Looks like we just need to remove "nyts" from the year variable that we created from the names of the tibbles in our list.
Here is our clean and wrangled data:
Rows: 95,465
Columns: 41
$ year <dbl> 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, …
$ psu <chr> "015438", "015438", "015438", "015438", "015438"…
$ finwgt <dbl> 216.7268, 324.9620, 324.9620, 397.1552, 264.8745…
$ stratum <chr> "BR3", "BR3", "BR3", "BR3", "BR3", "BR3", "BR3",…
$ Age <fct> 18, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, …
$ Sex <fct> female, male, male, male, female, female, male, …
$ Grade <fct> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, …
$ ECIGT <lgl> FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRU…
$ ECIGAR <lgl> TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, FA…
$ ESLT <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EELCIGT <lgl> FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRU…
$ EROLLCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EFLAVCIGTS <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, F…
$ EBIDIS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ EFLAVCIGAR <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, TRUE, F…
$ EHOOKAH <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ EPIPE <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ ESNUS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, …
$ EDISSOLV <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CCIGT <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, …
$ CCIGAR <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, …
$ CSLT <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CELCIGT <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, …
$ CROLLCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CFLAVCIGTS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CBIDIS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CHOOKAH <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CPIPE <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CSNUS <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ CDISSOLV <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,…
$ brand_ecig <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ menthol <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ clove_spice <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ fruit <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ chocolate <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ alcoholic_drink <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ candy_dessert_sweets <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ other <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ no_use <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ EHTP <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ CHTP <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …Reminder: Current users are a subset of ever users.
For lots of pivot_longer() that are the same… can use pivot_longer_spec()
spec <- relig_income %>% build_longer_spec( cols = -religion, names_to = “income”, values_to = “count” ) pivot_longer_spec(relig_income, spec)
But need to decide if we are going to create some new data frames…
plot1 <- nyts_data %>%
mutate(tobacco_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE),
tobacco_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(tobacco_ever = case_when(tobacco_sum_ever > 0 ~ TRUE,
tobacco_sum_ever ==0 ~ FALSE),
tobacco_current = case_when(tobacco_sum_current > 0 ~ TRUE,
tobacco_sum_current ==0 ~ FALSE)) %>%
group_by(year) %>%
summarise(tobacco_ever_year=(sum(tobacco_ever, na.rm = TRUE)*100)/
sum(!is.na(tobacco_ever)),
tobacco_current_year=(sum(tobacco_current, na.rm = TRUE)*100)/
sum(!is.na(tobacco_current))) %>%
rename("Ever"=tobacco_ever_year,
"Current"=tobacco_current_year) %>%
pivot_longer(cols = -year, names_to = "User", values_to = "Percentage of students")%>%
ggplot(aes(x=year,y=`Percentage of students`, linetype=User)) +
geom_line() +
geom_point(show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0,70,by=10),
labels = seq(0,70,by=10),
limits = c(0,70)) +
theme_minimal() +
theme(legend.position = "bottom",
axis.title.x = element_blank()) +
labs(title = "How does nicotine use vary over the years?",
subtitle = "Current and ever users of nicotine products",
y = "% of students")
plot1 
plot2 <- nyts_data %>%
group_by(year,
Sex) %>%
summarise(EELCIGT_year=(sum(EELCIGT, na.rm = TRUE)*100)/
sum(!is.na(EELCIGT)),
CELCIGT_year=(sum(CELCIGT, na.rm = TRUE)*100)/
sum(!is.na(CELCIGT))) %>%
filter(!is.na(Sex)) %>%
rename("E-cigarettes, Ever"=EELCIGT_year,
"E-cigarettes, Current"=CELCIGT_year) %>%
#converting all columns between and including `E-cigarettes, Ever` and `E-cigarettes, Current` into one column called category
pivot_longer(cols = `E-cigarettes, Ever`:`E-cigarettes, Current`, names_to = "Category", values_to = "Percentage of students")%>%
mutate(User = case_when(Category == "E-cigarettes, Ever" ~ "Ever",
Category == "E-cigarettes, Current" ~ "Current")) %>%
# mutate(Sex = case_when(female == TRUE ~ "Females",
# female == FALSE ~ "Males")) %>%
ggplot(aes(x=year,y=`Percentage of students`, color=Sex, linetype=User)) +
geom_line() +
geom_point(show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
theme_minimal() +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
axis.title.x = element_blank()) +
labs(title = "How do vaping rates compare between males and females?",
subtitle = "Current and ever users by gender",
y = "% of students")
plot2
What vaping brands and flavors appear to be used the most frequently?
Huang J, Duan Z, Kwok J, et al. Tob Control 2019;28:146–151.
plot3 <- nyts_data %>%
filter(year==2019) %>%
group_by(brand_ecig) %>%
filter(!is.na(brand_ecig)) %>%
summarise(n = n()) %>%
mutate(total = sum(n),
Percent = n*100/total) %>%
mutate(brand_ecig = fct_reorder(brand_ecig, desc(Percent))) %>%
ggplot(aes(x=brand_ecig,y=Percent, fill=brand_ecig)) +
geom_bar(stat="identity") +
theme_minimal() +
theme(legend.position = "none",
axis.text.x = element_text(size=10),
axis.title.x = element_blank()) +
labs(title = "What vaping brands appear to be used the most frequently?",
subtitle = "Brand of e-cigarette most frequently used in the last 30 days (2019)",
y = "% of e-cigarette users responding")
plot3
plot4 <- nyts_data %>%
filter(year!=2015) %>%
group_by(year) %>%
summarise(Menthol=(sum(menthol, na.rm = TRUE)*100)/
sum(!is.na(menthol)),
`Clove or Spice`=(sum(clove_spice, na.rm = TRUE)*100)/
sum(!is.na(clove_spice)),
`Fruit`=(sum(fruit, na.rm = TRUE)*100)/
sum(!is.na(fruit)),
`Chocolate`=(sum(chocolate, na.rm = TRUE)*100)/
sum(!is.na(chocolate)),
`Alcoholic Drink`=(sum(alcoholic_drink, na.rm = TRUE)*100)/
sum(!is.na(alcoholic_drink)),
`Candy/Desserts/Sweets`=(sum(candy_dessert_sweets, na.rm = TRUE)*100)/sum(!is.na(candy_dessert_sweets))) %>%
pivot_longer(cols = -year, names_to = "Percentage of students", values_to = "Flavor") %>%
ggplot(aes(x=year, y=`Percentage of students`, color=Flavor)) +
geom_line() +
geom_point(show.legend = FALSE) +
theme_minimal() +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
axis.title.x = element_blank()) +
labs(title = "What flavors appear to be used the most frequently in nicotine products?",
subtitle = "Flavors of tobacco products used in the past 30 days")
plot4 
plot5 <- nyts_data %>%
filter(year!=2015) %>%
filter(menthol==TRUE|
clove_spice==TRUE|
fruit==TRUE|
chocolate==TRUE|
alcoholic_drink==TRUE|
candy_dessert_sweets==TRUE|
other==TRUE) %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
mutate(ecig_only_ever = case_when(ecig_ever == TRUE &
non_ecig_ever ==FALSE ~ TRUE,
TRUE ~ FALSE),
ecig_only_current = case_when(ecig_current == TRUE &
non_ecig_ever ==FALSE ~ TRUE,
TRUE ~ FALSE),
non_ecig_only_ever = case_when(non_ecig_ever == TRUE &
ecig_ever ==FALSE ~ TRUE,
TRUE ~ FALSE),
non_ecig_only_current = case_when(non_ecig_current == TRUE &
ecig_ever ==FALSE ~ TRUE,
TRUE ~ FALSE)) %>%
mutate(Group = case_when(ecig_only_ever==TRUE |
ecig_only_current==TRUE ~ "Only e-cigarettes",
non_ecig_only_ever==TRUE |
non_ecig_only_current==TRUE ~ "Only other products",
TRUE ~ "Both")) %>%
filter(Group!="Both") %>%
group_by(year, Group) %>%
summarise(`Menthol`=(sum(menthol, na.rm = TRUE)*100)/
sum(!is.na(menthol)),
`Clove or Spice`=(sum(clove_spice, na.rm = TRUE)*100)/
sum(!is.na(clove_spice)),
`Fruit`=(sum(fruit, na.rm = TRUE)*100)/sum(!is.na(fruit)),
`Chocolate`=(sum(chocolate, na.rm = TRUE)*100)/
sum(!is.na(chocolate)),
`Alcoholic Drink`=(sum(alcoholic_drink, na.rm = TRUE)*100)/
sum(!is.na(alcoholic_drink)),
`Candy/Desserts/Sweets`=(sum(candy_dessert_sweets, na.rm = TRUE)*100)/
sum(!is.na(candy_dessert_sweets)),
`Other`=(sum(other, na.rm = TRUE)*100)/
sum(!is.na(other)),
Respondents=n()) %>%
#converting all columns between and including Menthol and Other to one column called Flavor
pivot_longer(cols = Menthol:Other, names_to = "Flavor", values_to = "Percentage of students") %>%
filter(!is.na(`Percentage of students`),
Flavor!="Other") %>%
group_by(Flavor) %>%
mutate(affirmative=(Respondents * `Percentage of students`)/100) %>%
mutate(flavor_mean = sum(affirmative)/sum(Respondents)) %>%
ungroup() %>%
mutate(flavor_mean_rank = dense_rank(flavor_mean),
Flavor = fct_reorder(Flavor, flavor_mean_rank)) %>%
ggplot(aes(x=year, y=`Percentage of students`, color=Group)) +
facet_wrap(.~Flavor,ncol=3) +
geom_line() +
geom_point(show.legend = FALSE) +
theme_minimal() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90)) +
labs(title = "Among users of only one type of product, what vaping flavors appear to be used the most frequently?",
subtitle = "Percent reporting only e-cigarette use vs only other nicotine product use")
plot5
Have vaping rates possibly influenced tobacco/nicotine use?
plot6 <- nyts_data %>%
group_by(year) %>%
summarise(ECIGT_year=(sum(ECIGT, na.rm = TRUE)*100)/
sum(!is.na(ECIGT)),
EELCIGT_year=(sum(EELCIGT, na.rm = TRUE)*100)/
sum(!is.na(EELCIGT)),
CCIGT_year=(sum(CCIGT, na.rm = TRUE)*100)/
sum(!is.na(CCIGT)),
CELCIGT_year=(sum(CELCIGT, na.rm = TRUE)*100)/
sum(!is.na(CELCIGT))) %>%
rename("Cigarettes, Ever"=ECIGT_year,
"E-cigarettes, Ever"=EELCIGT_year,
"Cigarettes, Current"=CCIGT_year,
"E-cigarettes, Current"=CELCIGT_year) %>%
pivot_longer(cols= -year, names_to = "Category", values_to = "Percentage of students")%>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(User = case_when(Category == "Cigarettes, Ever" ~ "Ever",
Category == "E-cigarettes, Ever" ~ "Ever",
Category == "Cigarettes, Current" ~ "Current",
Category == "E-cigarettes, Current" ~ "Current")) %>%
mutate(Product = case_when(Category == "Cigarettes, Ever" ~ "Cigarettes",
Category == "E-cigarettes, Ever" ~ "E-cigarettes",
Category == "Cigarettes, Current" ~ "Cigarettes",
Category == "E-cigarettes, Current" ~ "E-cigarettes")) %>%
ggplot(aes(x=year,y=`Percentage of students`, color=Product, linetype=User)) +
geom_line() +
geom_point(show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
theme_minimal() +
theme(legend.position = "bottom",
axis.title.x = element_blank()) +
labs(title = "How does e-cigarette use compare to cigarette use?",
subtitle = "Current and ever users of e-cigarettes and cigarettes",
y = "% of students")
plot6
plot7 <- nyts_data %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
group_by(year) %>%
summarise(ecig_ever_year=(sum(ecig_ever, na.rm = TRUE)*100)/
sum(!is.na(ecig_ever)),
ecig_current_year=(sum(ecig_current, na.rm = TRUE)*100)/
sum(!is.na(ecig_current)),
non_ecig_ever_year=(sum(non_ecig_ever, na.rm = TRUE)*100)/
sum(!is.na(non_ecig_ever)),
non_ecig_current_year=(sum(non_ecig_current, na.rm = TRUE)*100)/
sum(!is.na(non_ecig_current))) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students") %>%
mutate(User = case_when(Category =="ecig_ever_year" ~ "Ever",
Category =="non_ecig_ever_year" ~ "Ever",
Category =="ecig_current_year" ~ "Current",
Category =="non_ecig_current_year" ~ "Current")) %>%
mutate(Product = case_when(Category =="ecig_ever_year" ~ "E-cigarettes",
Category =="non_ecig_ever_year" ~ "Other products",
Category =="ecig_current_year" ~ "E-cigarettes",
Category =="non_ecig_current_year" ~ "Other products")) %>%
filter(User=="Ever") %>%
ggplot(aes(x=year,y=`Percentage of students`, color=Product)) +
geom_line(linetype=1) + # geom_bar(stat="identity", position = "dodge", color="black") +
geom_point(show.legend = FALSE) +
scale_y_continuous(breaks = seq(10, 60, by = 10), limits = c(10,60)) +
theme_minimal() +
theme(legend.position = "bottom",
axis.title.x = element_blank()) +
labs(title = "How does e-cigarette ever use compare to ever use of other products over the years?",
subtitle = "E-cigarette and non-e-cigarette use",
y = "% of students")
plot7
plot8 <- nyts_data %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
group_by(year) %>%
summarise(ecig_ever_year=(sum(ecig_ever, na.rm = TRUE)*100)/
sum(!is.na(ecig_ever)),
ecig_current_year=(sum(ecig_current, na.rm = TRUE)*100)/
sum(!is.na(ecig_current)),
non_ecig_ever_year=(sum(non_ecig_ever, na.rm = TRUE)*100)/
sum(!is.na(non_ecig_ever)),
non_ecig_current_year=(sum(non_ecig_current, na.rm = TRUE)*100)/
sum(!is.na(non_ecig_current))) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students") %>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(User = case_when(Category =="ecig_ever_year" ~ "Ever",
Category =="non_ecig_ever_year" ~ "Ever",
Category =="ecig_current_year" ~ "Current",
Category =="non_ecig_current_year" ~ "Current")) %>%
mutate(Product = case_when(Category =="ecig_ever_year" ~ "E-cigarettes",
Category =="non_ecig_ever_year" ~ "Other products",
Category =="ecig_current_year" ~ "E-cigarettes",
Category =="non_ecig_current_year" ~ "Other products")) %>%
ggplot(aes(x=year,y=`Percentage of students`, color=Product, linetype=User)) +
geom_line() +
geom_point(show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0, 60, by = 10), limits = c(0,60)) +
theme_minimal() +
theme(legend.position = "bottom",
axis.title.x = element_blank()) +
labs(title = "How does e-cigarette use compare to use of other products over the years?",
subtitle = "E-cigarette and non-e-cigarette use",
y = "% of students")
plot8
plotA_uw <- plot1 +
theme(axis.title.x = element_blank(),
legend.position = "none") +
labs(title = "Nicotine product users more prevalent after 2017",
subtitle = NULL,
y = "% of students")
plotB_uw <- plot7 +
theme(axis.title.x = element_blank(),
legend.position = "none") +
labs(title = "% Ever trying e-cigarettes increases &\n% ever trying other products decreases",
subtitle = NULL,
y = "% of students")
plotC_uw <- plot8 +
theme(axis.title.x = element_blank(),
legend.position = "none") +
labs(title = "% Using e-cigarettes increases &\n% using Other products decreases",
subtitle = NULL,
y = "% of students")
title_uw <- ggdraw() +
draw_label(
"Have vaping rates possibly influenced tobacco/nicotine use?",
fontface = 'bold',
size=14,
x = 0,
hjust = 0
) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsA_uw <- plot_grid(plotA_uw,
rel_widths = c(1,1))
plotsBC_uw <- plot_grid(plotB_uw,
plotC_uw,
rel_widths = c(1,1))
legend_uw <- get_legend(plotB_uw +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
figure_uw <- plot_grid(title_uw,
plotsA_uw,
plotsBC_uw,
legend_uw,
ncol = 1,
rel_heights = c(0.1,
1,
1,
0.1),
scale = 1.0)
figure_uw
plotA_w <- nyts_data %>%
mutate(tobacco_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE),
tobacco_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(tobacco_ever = case_when(tobacco_sum_ever > 0 ~ TRUE,
tobacco_sum_ever ==0 ~ FALSE),
tobacco_current = case_when(tobacco_sum_current > 0 ~ TRUE,
tobacco_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt, nest=TRUE) %>%
group_by(year) %>%
summarise(tobacco_ever_year = survey_mean(tobacco_ever, vartype = "ci", na.rm=TRUE),
tobacco_current_year = survey_mean(tobacco_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
pivot_longer(cols = -year, names_to = "Type", values_to = "Percentage of students") %>%
# gather(key=Type,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Type) ~ "Lower",
grepl("_upp", Type) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("ever", Type) ~ "Ever",
grepl("current", Type) ~ "Current",
TRUE ~ "Mean")) %>%
dplyr::select(-Type) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
ggplot(aes(x=year,y=Mean)) +
geom_line(aes(linetype=User)) +
geom_linerange(aes(ymin = Lower, ymax = Upper), show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0,70,by=10),
labels = seq(0,70,by=10),
limits = c(0,70)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "Nicotine product users more prevalent after 2017",
y = "% of students")
plotB_w <- nyts_data %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt, nest=TRUE) %>%
group_by(year) %>%
summarise(ecig_ever_year = survey_mean(ecig_ever, vartype = "ci", na.rm=TRUE),
ecig_current_year = survey_mean(ecig_current, vartype = "ci", na.rm=TRUE),
non_ecig_ever_year = survey_mean(non_ecig_ever, vartype = "ci", na.rm=TRUE),
non_ecig_current_year = survey_mean(non_ecig_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
dplyr::select(year,
ecig_ever_year,
ecig_current_year,
non_ecig_ever_year,
non_ecig_current_year,
contains("low"),
contains("upp")) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students") %>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Category) ~ "Lower",
grepl("_upp", Category) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("current", Category) ~ "Current",
TRUE ~ "Ever",),
Product = case_when(grepl("non_ecig", Category) ~ "Other products",
TRUE ~ "E-cigarettes")) %>%
dplyr::select(-Category) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
filter(User=="Ever") %>%
dplyr::rename("Lower_temp" = Upper,
"Upper_temp" = Lower) %>%
dplyr::rename("Lower"=Lower_temp,
"Upper"=Upper_temp) %>%
ggplot(aes(x=year,y=Mean, color=Product)) +
geom_line(linetype=1) +
geom_linerange(aes(ymin = Lower, ymax = Upper), show.legend = FALSE) +
scale_y_continuous(breaks = seq(10, 60, by = 10), limits = c(10,60)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "% ever trying e-cigarettes increases &\n% ever trying other products decreases",
y = "% of students")
#### the wrangling looks the same as the above plot...
plotC_w <- nyts_data %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt, nest=TRUE) %>%
group_by(year) %>%
summarise(ecig_ever_year = survey_mean(ecig_ever, vartype = "ci", na.rm=TRUE),
ecig_current_year = survey_mean(ecig_current, vartype = "ci", na.rm=TRUE),
non_ecig_ever_year = survey_mean(non_ecig_ever, vartype = "ci", na.rm=TRUE),
non_ecig_current_year = survey_mean(non_ecig_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
dplyr::select(year,
ecig_ever_year,
ecig_current_year,
non_ecig_ever_year,
non_ecig_current_year,
contains("low"),
contains("upp")) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students") %>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Category) ~ "Lower",
grepl("_upp", Category) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("ever", Category) ~ "Ever",
grepl("current", Category) ~ "Current"),
Product = case_when(grepl("non_ecig", Category) ~ "Other products",
TRUE ~ "E-cigarettes")) %>%
dplyr::select(-Category) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
ggplot(aes(x=year,y=Mean, color=Product)) +
geom_line(aes(linetype=User)) +
geom_linerange(aes(ymin = Lower, ymax = Upper), show.legend = FALSE) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0, 60, by = 10), limits = c(0,60)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "% Using e-cigarettes increases &\n% using Other products decreases",
y = "% of students")
title_w <- ggdraw() +
draw_label(
expression("Have vaping rates possibly influenced tobacco/nicotine use?"),
fontface = 'bold',
size=14,
x = 0,
hjust = 0
) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsA_w <- plot_grid(plotA_w,
rel_widths = c(1),
align = "v",
axis = "bt")
plotsBC_w <- plot_grid(plotB_w,
plotC_w,
rel_widths = c(1,1),
align = "v",
axis = "bt")
legend_w <- get_legend(plotB_w +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
figure_w <- plot_grid(title_w,
plotsA_w,
plotsBC_w,
legend_w,
ncol = 1,
rel_heights = c(0.1,
1,
1,
0.1),
scale = 1.0)
figure_w
plotA_w_8 <- nyts_data %>%
filter((Grade == "8" & year == 2015) |
(Grade == "9" & year == 2016) |
(Grade == "10" & year == 2017) |
(Grade == "11" & year == 2018) |
(Grade == "12" & year == 2019)
) %>%
mutate(tobacco_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE),
tobacco_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(tobacco_ever = case_when(tobacco_sum_ever > 0 ~ TRUE,
tobacco_sum_ever ==0 ~ FALSE),
tobacco_current = case_when(tobacco_sum_current > 0 ~ TRUE,
tobacco_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt) %>%
group_by(year) %>%
summarise(tobacco_ever_year = survey_mean(tobacco_ever, vartype = "ci", na.rm=TRUE),
tobacco_current_year = survey_mean(tobacco_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
pivot_longer(cols = -year, names_to = "Type", values_to = "Percentage of students")%>%
# gather(key=Type,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Type) ~ "Lower",
grepl("_upp", Type) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("ever", Type) ~ "Ever",
grepl("current", Type) ~ "Current",
TRUE ~ "Mean")) %>%
dplyr::select(-Type) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
ggplot(aes(x=year,y=Mean)) +
geom_line(aes(linetype=User)) +
geom_linerange(aes(ymin = Lower, ymax = Upper)) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0,70,by=10),
labels = seq(0,70,by=10),
limits = c(0,70)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "Nicotine product users becoming increasingly prevalent",
y = "% of students")
plotB_w_8 <- nyts_data %>%
filter((Grade == "8" & year == 2015) |
(Grade == "9" & year == 2016) |
(Grade == "10" & year == 2017) |
(Grade == "11" & year == 2018) |
(Grade == "12" & year == 2019)
) %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt) %>%
group_by(year) %>%
summarise(ecig_ever_year = survey_mean(ecig_ever, vartype = "ci", na.rm=TRUE),
ecig_current_year = survey_mean(ecig_current, vartype = "ci", na.rm=TRUE),
non_ecig_ever_year = survey_mean(non_ecig_ever, vartype = "ci", na.rm=TRUE),
non_ecig_current_year = survey_mean(non_ecig_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
dplyr::select(year,
ecig_ever_year,
ecig_current_year,
non_ecig_ever_year,
non_ecig_current_year,
contains("low"),
contains("upp")) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students")%>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Category) ~ "Lower",
grepl("_upp", Category) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("current", Category) ~ "Current",
TRUE ~ "Ever",),
Product = case_when(grepl("non_ecig", Category) ~ "Other products",
TRUE ~ "E-cigarettes")) %>%
dplyr::select(-Category) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
filter(User=="Ever") %>%
dplyr::rename("Lower_temp" = Upper,
"Upper_temp" = Lower) %>%
dplyr::rename("Lower"=Lower_temp,
"Upper"=Upper_temp) %>%
ggplot(aes(x=year,y=Mean, color=Product)) +
geom_line(linetype=1) +
geom_linerange(aes(ymin = Lower, ymax = Upper)) +
scale_y_continuous(breaks = seq(10, 60, by = 10), limits = c(10,60)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "% ever trying nicotine products increases",
y = "% of students")
plotC_w_8 <- nyts_data %>%
filter((Grade == "8" & year == 2015) |
(Grade == "9" & year == 2016) |
(Grade == "10" & year == 2017) |
(Grade == "11" & year == 2018) |
(Grade == "12" & year == 2019)
) %>%
mutate(ecig_sum_ever = select(., EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
ecig_sum_current = select(., CELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_ever = select(., starts_with("E", ignore.case = FALSE)) %>%
select(.,-EELCIGT) %>%
apply(1, sum, na.rm=TRUE),
non_ecig_sum_current = select(., starts_with("C", ignore.case = FALSE)) %>%
select(.,-CELCIGT) %>%
apply(1, sum, na.rm=TRUE)) %>%
mutate(ecig_ever = case_when(ecig_sum_ever > 0 ~ TRUE,
ecig_sum_ever ==0 ~ FALSE),
ecig_current = case_when(ecig_sum_current > 0 ~ TRUE,
ecig_sum_current ==0 ~ FALSE),
non_ecig_ever = case_when(non_ecig_sum_ever > 0 ~ TRUE,
non_ecig_sum_ever ==0 ~ FALSE),
non_ecig_current = case_when(non_ecig_sum_current > 0 ~ TRUE,
non_ecig_sum_current ==0 ~ FALSE)) %>%
as_survey_design(strata = stratum, ids = psu, weight = finwgt) %>%
group_by(year) %>%
summarise(ecig_ever_year = survey_mean(ecig_ever, vartype = "ci", na.rm=TRUE),
ecig_current_year = survey_mean(ecig_current, vartype = "ci", na.rm=TRUE),
non_ecig_ever_year = survey_mean(non_ecig_ever, vartype = "ci", na.rm=TRUE),
non_ecig_current_year = survey_mean(non_ecig_current, vartype = "ci", na.rm=TRUE)) %>%
mutate_at(vars(-year), "*", 100) %>%
dplyr::select(year,
ecig_ever_year,
ecig_current_year,
non_ecig_ever_year,
non_ecig_current_year,
contains("low"),
contains("upp")) %>%
pivot_longer(cols = -year, names_to = "Category", values_to = "Percentage of students")%>%
# gather(key=Category,
# value=`Percentage of students`,
# -year) %>%
mutate(Estimate = case_when(grepl("_low", Category) ~ "Lower",
grepl("_upp", Category) ~ "Upper",
TRUE ~ "Mean"),
User = case_when(grepl("ever", Category) ~ "Ever",
grepl("current", Category) ~ "Current"),
Product = case_when(grepl("non_ecig", Category) ~ "Other products",
TRUE ~ "E-cigarettes")) %>%
dplyr::select(-Category) %>%
pivot_wider(names_from = Estimate,
values_from = `Percentage of students`) %>%
#spread(Estimate, `Percentage of students`) %>%
ggplot(aes(x=year,y=Mean, color=Product)) +
geom_line(aes(linetype=User)) +
geom_linerange(aes(ymin = Lower, ymax = Upper)) +
scale_linetype_manual(values = c(2,1)) +
scale_y_continuous(breaks = seq(0, 60, by = 10), limits = c(0,60)) +
theme_minimal() +
theme(legend.position = "none",
axis.title.x = element_blank()) +
labs(title = "E-cigarette use surpasses use of other nicotine products",
y = "% of students")
title_w_8 <- ggdraw() +
draw_label(
expression("Among"~8^th~"graders in 2015, have vaping rates possibly influenced tobacco/nicotine use?"),
fontface = 'bold',
size=14,
x = 0,
hjust = 0
) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsA_w_8 <- plot_grid(plotA_w_8,
rel_widths = c(1),
align = "v",
axis = "bt")
plotsBC_w_8 <- plot_grid(plotB_w_8,
plotC_w_8,
rel_widths = c(1,1),
axis = "bt")
legend_w_8 <- get_legend(plotB_w_8 +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
figure_w_8 <- plot_grid(title_w_8,
plotsA_w_8,
plotsBC_w_8,
legend_w_8,
ncol = 1,
rel_heights = c(0.1,
1,
1,
0.1),
scale = 1.0
)
figure_w_8
title_final <- ggdraw() +
draw_label(
expression("Have vaping rates possibly influenced tobacco/nicotine use?"),
fontface = 'bold',
size=16,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
subtitle_uw_final <- ggdraw() +
draw_label(
expression(~6^th~"-"~12^th~"graders, unweighted"),
size=12,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
subtitle_w_final <- ggdraw() +
draw_label(
expression(~6^th~"-"~12^th~"graders, weighted"),
fontface = 'bold',
size=12,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
subtitle_w_8_final <- ggdraw() +
draw_label(
expression(~8^th~"graders in 2015, weighted"),
fontface = 'bold',
size=12,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
subtitle_final <- plot_grid(subtitle_uw_final,
subtitle_w_final,
subtitle_w_8_final,
ncol = 3)
plotsA_title_final <- ggdraw() +
draw_label(
expression("Prevalence of e-cigarette use by user type"),
size=14,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsA_final <- plot_grid(plotA_uw + theme(plot.title = element_blank()),
plotA_w + theme(plot.title = element_blank()),
plotA_w_8 + theme(plot.title = element_blank()),
ncol = 3,
align = "v",
axis = "bt")
plotsB_title_final <- ggdraw() +
draw_label(
expression("Prevalence of ever use by product type"),
size=14,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsB_final <- plot_grid(plotB_uw + theme(plot.title = element_blank()),
plotB_w + theme(plot.title = element_blank()),
plotB_w_8 + theme(plot.title = element_blank()),
ncol = 3,
align = "v",
axis = "bt")
plotsC_title_final <- ggdraw() +
draw_label(
expression("Prevalence of nicotine product use by product & user type"),
size=14,
x = 0.5) +
theme(
plot.margin = margin(0, 0, 0, 0)
)
plotsC_final <- plot_grid(plotC_uw + theme(plot.title = element_blank()),
plotC_w + theme(plot.title = element_blank()),
plotC_w_8 + theme(plot.title = element_blank()),
ncol = 3,
align = "v",
axis = "bt")
legend_final <- get_legend(plotB_w +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
final_plot <- plot_grid(title_final,
plotsA_title_final,
subtitle_final,
plotsA_final,
plotsB_title_final,
subtitle_final,
plotsB_final,
plotsC_title_final,
subtitle_final,
plotsC_final,
legend_final,
ncol = 1,
rel_heights = c(0.2,
0.2,
0.1,
1,
0.2,
0.1,
1,
0.2,
0.1,
1,
0.1))
final_plot
final_plot above for a different cohort of your choice.Ever and current variables are limited to those shared by all years of data included in this case study.
New code: Knit time: 32.742 secs. Previous code: ~ 3 m.
I had difficulty producing a plot that succinctly presented a trend. It’s very easy to produce plots that are very useful once one is familiar with the data. Some plots, however, cannot stand alone and need additional context to be clear for those without prior knowledge of the data. When I first shared a plot I had been working on with others, it became clear that in my effort to present a complicated idea briefly I had left out information that would make the trend easily interpretable. To solve this issue, I began to present visualizations of the data alongside my original plot. The final figure I created contained several additional plots, each presenting the same trend at a different level than my initial plot.
My “centerpiece” plot is the middle plot in final_plot. The 8 plots around it help provide a very clear picture of what is going on in the US with regards to e-cigarette use and nicotine product use at large. On its own, it’s difficult to understand the trends in the US and how important the weighting scheme is for inference. Once you add the left and right columns, it’s clear what is going on.
R version 4.0.1 (2020-06-06)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.5
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] cowplot_1.0.0 srvyr_0.3.10 forcats_0.5.0 ggplot2_3.3.1 tidyr_1.1.0
[6] readr_1.3.1 dplyr_1.0.0 purrr_0.3.4 stringr_1.4.0 readxl_1.3.1
[11] knitr_1.28 here_0.1
loaded via a namespace (and not attached):
[1] Rcpp_1.0.4.6 highr_0.8 cellranger_1.1.0 pillar_1.4.4
[5] compiler_4.0.1 base64enc_0.1-3 tools_4.0.1 digest_0.6.25
[9] lattice_0.20-41 evaluate_0.14 lifecycle_0.2.0 tibble_3.0.1
[13] gtable_0.3.0 pkgconfig_2.0.3 rlang_0.4.6 Matrix_1.2-18
[17] cli_2.0.2 DBI_1.1.0 yaml_2.2.1 xfun_0.14
[21] withr_2.2.0 mitools_2.4 generics_0.0.2 vctrs_0.3.0
[25] hms_0.5.3 rprojroot_1.3-2 grid_4.0.1 tidyselect_1.1.0
[29] glue_1.4.1 R6_2.4.1 fansi_0.4.1 survival_3.1-12
[33] rmarkdown_2.2 farver_2.0.3 magrittr_1.5 splines_4.0.1
[37] backports_1.1.7 scales_1.1.1 ellipsis_0.3.1 htmltools_0.4.0
[41] assertthat_0.2.1 colorspace_1.4-1 labeling_0.3 utf8_1.1.4
[45] stringi_1.4.6 survey_4.0 munsell_0.5.0 crayon_1.3.4